UNCLASSIFIED 


Defense  Technical  Information  Center 
Compilation  Part  Notice 

ADP012017 

TITLE:  Curves  from  Motion,  Motion  from  Curves 
DISTRIBUTION:  Approved  for  public  release,  distribution  unlimited 


This  paper  is  part  of  the  following  report: 

TITLE:  International  Conference  on  Curves  and  Surfaces  [4th],  Saint-Malo, 
France,  1-7  July  1999.  Proceedings,  Volume  1.  Curve  and  Surface  Design 

To  order  the  complete  compilation  report,  use:  ADA399461 

The  component  part  is  provided  here  to  allow  users  access  to  individually  authored  sections 
of  proceedings,  annals,  symposia,  etc.  However,  the  component  should  be  considered  within 
the  context  of  the  overall  compilation  report  and  not  as  a stand-alone  technical  report. 

The  following  component  part  numbers  comprise  the  compilation  report: 

ADPO 1 20 1 0 thru  ADPO 1 2054 


UNCLASSIFIED 


Curves  from  Motion,  Motion  from  Curves 


Rida  T.  Farouki 


Abstract.  Geometry  and  kinematics  have  been  intimately  connected  in 
their  historical  evolution  and,  although  it  is  currently  less  fashionable,  the 
further  development  of  such  connections  is  crucial  to  many  computer-aided 
design  and  manufacturing  applications.  In  this  survey,  we  explore  a variety 
of  classical  and  modern  problems  that  illustrate  how  simple  rules  of  motion 
produce  interesting  curves  and,  conversely,  the  computational  problems  of 
generating  motions  with  prescribed  paths  and  speeds.  These  encompass 
the  geometry  of  trajectories  under  centripetal  forces;  the  transformation  of 
rotary  motion  into  motion  along  general  curves  by  mechanisms;  real-time 
curve  interpolators  for  digital  motion  control;  and  the  description  of  spatial 
motions  that  involve  variations  of  both  position  and  orientation.  Such  case 
studies  illustrate  some  of  the  intellectual  appeal,  and  practical  importance, 
of  a sustained  dialog  between  the  study  of  curves  and  of  motions. 


§1.  Preamble 

Our  intent  in  this  paper  is  to  survey  the  intricate  web  of  historical  connections 
between  geometry  and  kinematics,  a theme  that  has  played  a key  role  in  the 
development  of  mechanics  and  analysis.  In  contemplating  this  theme,  we  are 
obviously  confronted  by  a profusion  of  interesting  and  fruitful  topics  — and 
we  are  thus  obliged  to  adopt  a rather  anecdotal  approach. 

Apart  from  its  intrinsic  interest,  we  choose  this  subject  with  the  hope  of 
promoting  greater  synergy  between  modern-day  problems  of  geometric  design 
and  motion  control.  Modern  CAD  systems  are  mainly  concerned  with  creating 
“static”  geometrical  descriptions  of  artifacts,  but  the  processes  by  which  these 
artifacts  are  actually  fabricated  often  involve  complicated  motions  of  a tool  — 
e.g.,  a cutter  in  a milling  machine,  or  a wire  electrode  in  electrical  discharge 
machining  — relative  to  a workpiece.  Compared  to  the  sophistication  of  CAD 
models,  current  methods  for  motion  planning  in  manufacturing  processes  are 
often  rather  crude  and  naive.  Thus,  there  is  much  scope  for  securing  greater 
precision  and  reliability,  through  the  use  of  advanced  mathematical  methods, 
in  the  relatively  undeveloped  field  of  manufacturing  geometry. 

The  symbiosis  between  geometry  and  kinematics  has  deep  historical  roots. 
Newton,  in  his  Quadrature  of  Curves  (1676),  aptly  characterizes  it  as  follows: 
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Lines  (curves)  are  described,  and  thereby  generated,  not  by  the 
apposition  of  parts  but  by  the  continued  motion  of  points  . . . These 
geneses  really  take  place  in  the  nature  of  things,  and  are  daily  seen 
in  the  motion  of  bodies. 

However,  this  has  not  always  been  a happy  union.  Insofar  as  it  embodies  both 
spatial  and  temporal  information,  kinematics  subsumes  geometry.  To  upgrade 
a curve  into  a motion  requires  the  ability  to  rectify  (or  measure  arc  lengths  of) 
curves.  As  a basic  philosophical  tenet,  Descartes  held  this  to  be  impossible  — 
see  §3  below  — and  he  sought  to  banish  all  curves  whose  definitions  explicitly 
or  implicitly  assume  rectifications  from  the  “rigorous”  domain  of  geometry  to 
the  nascent  (and  less-exact)  science  of  mechanics.  Although,  in  modern  times, 
the  philosophical/existential  problem  of  rectification  is  no  longer  troublesome, 
we  must  still  address  the  computational  difficulties  it  entails  (see  §3). 

The  antithesis  of  Descartes’  attempt  to  divorce  geometry  from  kinematics 
would  ultimately  find  its  logical  expression,  in  the  context  of  the  special  theory 
of  relativity,  with  Minkowski’s  introduction  [42]  of  the  concept  of  “space-time” 
as  the  most  natural  setting  for  the  study  of  physical  phenomena: 

Henceforth  space  by  itself,  and  time  by  itself,  are  doomed  to 
fade  away  into  mere  shadows,  and  only  a kind  of  union  of  the  two 
will  preserve  an  independent  reality. 

The  recent  introduction  of  the  Minkowski  metric  of  space-time  into  problems 
of  geometric  design  reveals  a remarkable  confluence  of  ideas  concerning  medial 
axis  transforms,  Pythagorean  hodographs,  envelopes,  and  offset  curves  [43], 
In  this  survey  we  shall  attempt,  through  a series  of  anecdotal  sketches,  to 
promote  greater  interest  in  the  relationship  between  geometry  and  kinematics, 
and  its  application  to  CAD/CAM  problems.  We  commence  in  §2  and  §3  with 
a brief  review  of  the  manner  in  which  curves  may  be  defined,  and  the  problems 
that  measurement  of  arc  length  incurs.  Perhaps  the  simplest  motion  is  that  of 
a particle  experiencing  a force  toward  a fixed  center,  of  magnitude  proportional 
to  a power  of  the  radial  distance  r.  As  is  well-known,  Newton  showed  that  an 
r~2  force  of  gravity  explains  the  conic  form  of  planetary  and  cometary  orbits. 
We  shall  see  in  §4,  however,  that  this  is  just  one  aspect  of  a more  profound 
theory  of  motion  under  centripetal  forces  in  Newton’s  Principia. 

Mechanisms  such  as  gears,  cams,  and  linkages  are  used  to  transform  forces 
and  motions  in  machinery.  In  §5  we  discuss  the  four-bar  linkage,  a mechanism 
that  directly  transforms  rotary  motion  into  motion  on  a general  curved  path. 
CNC  machines  offer  a more  flexible  approach  to  motion  generation,  based  on 
sophisticated  servo-systems  that  drive  linear  or  rotary  axes  in  a coordinated 
manner.  In  §6  we  discuss  the  problem  of  real-time  interpolators,  which  must 
accurately  and  efficiently  interpret  the  path  and  speed  information  to  generate 
“reference  point”  data  required  by  the  digital  control  algorithm. 

The  preceding  examples  are  concerned  with  motion  in  Euclidean  spaces. 
A motion  that  involves  not  just  positional  but  also  orientational  coordinates 
(such  as  in  5-axis  machining)  may  be  regarded  as  the  motion  of  a point  in  a 
higher-dimensional,  non-Euclidean  (soma)  space.  Some  subtle  problems  that 
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arise  with  such  motions  are  discussed  in  §7.  Finally,  §8  offers  some  concluding 
thoughts  on  our  theme  of  connections  between  geometry  and  kinematics. 

§2.  Curves  and  Motions 

Analytic  geometry  has  its  origins  in  the  computational  investigation  of  curves 
specified  by  suitable  coordinate  equations.  There  are  basically  two  ways  to 
define  a plane  curve  in  terms  of  Cartesian  coordinates  ( x , y).  We  may  select  a 
predicate  function,  that  indicates  whether  or  not  each  point  in  the  plane  lies 
on  the  curve  — this  is  typically  a bivariate  polynomial  / in  the  coordinates, 
and  the  curve  is  the  locus  of  points  on  which  the  polynomial  vanishes: 

f(x,y)  = 0.  (1) 

On  the  other  hand,  we  may  choose  a pair  of  generating  functions 

x(t) , y(t)  (2) 

that  produce  an  ordered  sequence  of  curve  points  when  evaluated  at  successive 
values  of  a continuous  “auxiliary  variable”  or  parameter  t. 

Whereas  the  implicit  description  (1)  is  essentially  “static,”  the  parametric 
form  (2)  offers  a more  “dynamic”  characterization  of  curves  — it  embodies  the 
suggestion  of  motion  along  a curve,  incurred  by  steady  increase  of  the  curve 
parameter  t.  It  is  a mistake,  however,  to  invest  too  much  hope  in  the  capacity 
of  parametric  curves  to  adequately  describe  motions.  Motion  specification  is 
concerned  as  much  with  the  instants  in  time  at  which  a body  assumes  given 
positions  along  a path,  and  corresponding  velocities  and  accelerations,  as  with 
the  path  geometry.  A motion  is  really  a geometrical  locus  in  Minkowski  space, 
with  one  temporal  and  one  or  more  spatial  dimensions. 

Of  course,  we  can  always  interpret  the  parameter  t as  time,  and  equations 
(2)  then  completely  specify  a motion.  However,  if  we  wish  to  use  only  “simple” 
(polynomial  or  rational)  functions,  such  motions  are  mathematical  curiosities: 
except  in  trivial  cases,  they  are  neither  solutions  to  appropriate  equations  of 
motion,  nor  do  they  represent  motions  of  practical  interest  that  we  may  wish 
to  impose  on  a given  locus.  To  emphasize  that  the  curve  parameter  generally 
lacks  any  geometrical  or  temporal  significance,  we  henceforth  use  the  Greek 
character  £ to  denote  it,  and  we  explicitly  reserve  t for  time. 

Connections  between  the  study  of  curves  and  motions  is  a recurrent  theme 
in  the  history  of  science  and  technology.  At  the  inception  of  analytic  geometry, 
motions  offered  an  intuitive  means  to  construct  and  analyze  loci  of  increasing 
sophistication:  see,  for  example,  the  remarkably  diverse  historical  applications 
of  the  roulettes  generated  by  the  rolling  motions  of  lines  and  circles  (cycloids, 
circle  involutes,  epicycloids  and  hypocycloids,  epitrochoids  and  hypotrochoids) 
discussed  in  [21].  Conversely,  the  conics  of  the  ancient  Greeks  make  a rather 
surprising  appearance  in  the  solution  to  the  premier  problem  of  dynamics:  the 
determination  of  planetary  orbits.  In  the  modern  computer  era,  the  problem 
of  producing  a desired  motion  along  a given  path  is  central  to  real-time  control 
of  manufacturing,  inspection,  robotic,  and  other  devices. 
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§3.  Towards  an  Impossible  Ideal 

Since  speed  on  a curved  path  is  the  rate  of  change  of  distance  with  time,  the 
problem  of  rectification  — i.e.,  the  measurement  of  arc  length  — is  evidently 
critical  to  the  description  of  motion.  This  problem,  however,  has  been  fraught 
with  computational  difficulty  since  Descartes  founded  analytic  geometry  in  an 
appendix  La  geometric  [15]  to  the  Discours  de  la  methode  pour  bien  conduire 
sa  raison  et  chercher  la  verite  dans  les  sciences  (1637).  He  asserts  that: 

Geometry  should  not  include  lines  (curves)  that  are  like  strings, 
in  that  they  are  sometimes  straight  and  sometimes  curved,  since  the 
ratios  between  straight  and  curved  lines  are  not  known,  and  I believe 
cannot  be  discovered  by  human  minds,  and  therefore  no  conclusion 
based  upon  such  ratios  can  be  accepted  as  rigorous  and  exact. 

Nevertheless,  Descartes’  dictum  began  to  crumble  almost  immediately  after 
its  enunciation,  amid  a flurry  of  counter-examples. 

For  example,  Galileo  [26]  realized  that,  when  a body  is  dropped  into  a hole 
drilled  through  the  center  of  a static  Earth,  it  executes  linear  simple  harmonic 
motion  across  the  full  Earth  diameter  under  the  influence  of  gravity: 

...  if  the  terrestrial  globe  were  perforated  through  the  center,  a 
cannon  ball  descending  through  the  hole  would  have  acquired  at  the 
center  such  an  impetus  from  its  speed  that  it  would  pass  beyond  the 
center  and  be  driven  upward  through  as  much  space  as  it  had  fallen, 
its  velocity  beyond  the  center  always  diminishing  with  losses  equal 
to  the  increments  acquired  in  the  descent  . . . 

On  a rotating  Earth,  however,  the  body  will  have  an  initial  tangential  velocity, 
and  the  nature  of  its  motion  in  the  hypothetical  case  of  “permeable”  matter 
(which  exerts  gravitational  forces  but  does  not  impede  motion)  is  not  obvious. 
Galileo’s  pupil,  Evangelista  Torricelli  (1608-1647),  conjectured  that  the  path 
would  be  a logarithmic  spiral  about  the  Earth’s  center,  described  by 

r = aeke  (3) 

in  polar  coordinates  (also  known  as  an  “equi-angular”  spiral,  since  the  tangent 
makes  a fixed  angle,  cot-1  k,  with  the  radius  vector).  Isaac  Newton  re-iterated 
this  conjecture  [2,60]  in  a letter  dated  November  28,  1679  to  Robert  Hooke, 
who  criticized  it  during  a Royal  Society  meeting  the  following  December  11. 
As  we  shall  see  in  §4  below,  Torricelli  and  Newton  were  quite  wrong:  the  path 
is  actually  — as  intuitively  argued  by  Hooke  — an  ellipse. 

During  his  investigations,  however,  Torricelli  discovered  a rectification  of 
the  spiral  (3)  in  1645  through  the  Archimedean  “method  of  exhaustion”  — he 
showed  that,  for  — oo  < 0 < 0,  the  arc  length  equals  the  length  of  the  tangent 
at  6 = 0 extended  to  the  y-axis  [9]  — namely,  %/l  + fc-2  a (see  Figure  1).  This 
is  a truly  remarkable  result,  since  the  curve  must  execute  an  infinite  number 
of  gyrations  about  the  origin  before  terminating  there! 
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Fig.  1.  The  arc  length  of  (3)  for  6 < 0 equals  the  length  PQ  on  the  tangent  line. 

The  logarithmic  spiral  was  also  known  to  Thomas  Harriot  (1560-1621)  as 
the  projection  of  a rhumb  line  on  the  Earth’s  surface  — i.e.,  the  path  traced  by 
a ship  that  sails  in  a fixed  compass  direction  — onto  the  equatorial  plane  [55], 
Jakob  Bernoulli  was  so  fascinated  by  the  self-similarity  of  this  curve  under 
coordinated  rotations  and  dilatations  about  the  origin,  that  he  arranged  to 
have  it  engraved  on  his  tombstone  with  the  caption  Eadem  mutata  resurgo  [3] 
— “Though  changed  I shall  arise  the  same.” 

Subsequently,  another  curve  was  rectified  by  Gilles  Personne  de  Roberval 
(1602-1675)  and  Christopher  Wren  (1632-1723)  — namely,  the  cycloid 

x{9)  = a{Q  — sm6),  y(9)  = o(l  — cos0)  (4) 

traced  by  a fixed  point  on  a circle  of  radius  a that  rolls  without  slipping  on  a 
straight  line  (see  Figure  2).  They  showed  that  a single  “arch”  (0  < 9 < 27r)  of 
this  curve  has  length  8a.  Although  it  has  now  fallen  into  obscurity,  the  cycloid 
was  a virtual  “proving  ground”  for  novel  mathematical  ideas  and  methods  in 
the  mid-17th  century:  it  caught  the  attention  of  all  the  leading  scientists,  and 
prompted  international  competitions  and  acrimonious  controversies.  See  [21] 
for  a discussion  of  its  tautochrone  and  brachistochrone  properties. 

To  Descartes,  however,  the  rectification  of  curves  such  as  the  spiral  (3)  and 
the  cycloid  (4)  was  suspect  — they  are  not  true  “geometrical”  (i.e.,  algebraic) 
but  rather  “mechanical”  (i.e.,  transcendental ) curves.  By  introducing  angular 
variables,  their  definitions  essentially  presuppose  a rectification  (of  the  circle). 
Nonetheless,  it  was  not  long  before  even  an  algebraic  curve,  under  the  scrutiny 
of  William  Neil  (1637-1670),  Hendrick  van  Heuraet  (1633-1660),  and  Pierre 
de  Fermat  (1601-1665),  succumbed  to  rectification  — the  cuspidal  cubic 

z(£)  = £2>  y(0  = (5) 

known  as  the  “semicubical  parabola.”  Its  arc  length  s,  measured  from  £ = 0, 
is  an  algebraic  function  of  the  parameter: 

m (9fc2£2  + 4)3/2  - 8 
S(e)  = 27P ' 
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Fig.  2.  Cycloid:  the  length  of  an  arch  is  eight  times  the  radius  of  the  rolling  circle. 

Ironically,  van  Heuraet  — an  associate  of  Huygens  — published  his  results  in 
an  appendix  to  van  Schooten’s  1659  Latin  version  of  Descartes,  Geometria  a 
Renato  Des  Cartes.  Neil’s  results  also  appeared  in  1659,  in  the  Tractatus  duo, 
prior  de  cycloide,  posterior  de  cissoide  published  by  John  Wallis,  and  Fermat’s 
work  followed  in  1660  in  De  linearum  curvarum  cum  lineis  rectis  comparatione 
dissertatio  geometrica  — an  appendix  to  a treatise  by  de  Lalouvere  (this  was 
the  only  publication  by  Fermat  to  appear  during  his  lifetime). 

Christiaan  Huygens  (1629-1695),  in  his  Horologium  oscillatorium of  1673, 
gave  a historical  account  [34]  of  these  rectifications  that  provoked  arguments 
over  the  priority  he  attributed  to  van  Heuraet  and  Wren  for  their  discovery  — 
see  Chapter  8 of  [31].  This  dispute  reflects  the  philosophical  importance  of  the 
rectification  problem,  which  had  been  held  impossible  through  long  tradition 
that  originated  with  Aristotle,  was  reinforced  in  the  11th  century  by  Ibn  Rushd 
(Averroes),  and  culminated  in  Descartes’  dogmatic  assertion.  Huygens’  theory 
of  evolutes  and  involutes,  employed  in  his  design  of  an  isochronous  pendulum 
clock  [21],  offered  profound  new  insight  into  this  age-old  problem.  The  cubic 
(5)  was  recognized  as  the  evolute  (locus  of  centers  of  curvature)  of  a parabola, 
while  the  cycloid  (4)  has  an  identical  (displaced)  cycloid  as  its  evolute. 

All  these  results  preceded  a formal  development  of  the  calculus.  Whereas 
the  latter  resolved  existential  issues  concerning  arc  lengths  by  defining  them, 
for  a (sufficiently  smooth)  parametric  curve  (x(£),  2/(0),  through  the  integral 

s{()  = f y/x'2(T)  + y'2(T)  dr  , (6) 

Jo 

there  remained  the  awkward  fact  that  this  does  not  admit  analytic  reduction 
except  in  trivial  or  exceptional  cases,  such  as  the  cubic  (5). 

In  fact,  with  the  emergence  of  differential  geometry,  it  became  customary 
to  assume  s = ( — i.e. , the  integrand  in  (6)  is  precisely  unity  — although  this 
natural  or  arc-length  parameterization  has  only  a hypothetical  existence:  it  is 
fundamentally  incompatible  with  curves  (except  straight  lines)  parameterized 
by  “simple”  functions.  This  fact  is  obvious  for  polynomial  curves,  but  its  proof 
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for  rational  curves  is  subtle  [23],  involving  Pythagorean  triples  of  polynomials, 
partial  fraction  decompositions,  and  the  calculus  of  residues. 

An  offshoot  to  this  proof  was  the  introduction  of  Pythagorean  hodograph 
(PH)  curves,  whose  hodograph  components  satisfy  the  condition 

AO  + 2/'2(0  = *2(0 

for  some  polynomial  <r(£),  and  are  thus  [22]  of  the  form 

*'(0  = AO  - v2(0 , 2 m = 2 f.(0«(0 , = AO  + v\0 

where  u(£),  u(£)  are  relatively  prime  polynomials.  For  PH  curves,  the  integral 
(6)  evidently  reduces  to  a polynomial  function  of  the  parameter  £.  This  fact 
proves  especially  propitious  in  the  formulation  of  real-time  CNC  interpolators 
for  digital  motion  control  applications  (see  §6  below). 

Venturing  beyond  PH  curves,  one  may  seek  to  encompass  a broader  class 
of  loci  by  allowing  more  complicated  arc-length  functions.  Suppose  we  allow 
s to  be  an  algebraic  function  of  the  parameter  f — i.e.,  there  exists  a bivariate 
polynomial  F( •,  ■)  such  that  (6)  satisfies 

f(*(0,0  = 0. 

An  algebraic  function  cannot,  in  general,  be  described  by  a simple  closed-form 
expression.  Nevertheless,  one  can  show  [53]  that  (6)  is  algebraic  if  and  only  if 
there  exists  a polynomial  h(£)  such  that 

[x,2(0  + y'2(0)h(0  = 

As  an  immediate  consequence,  if  the  function  (6)  is  algebraic,  it  must  have  the 
simple  form  s(£)  = 2 y/h(t;)  + constant  (note  that  the  PH  curves  are  subsumed 
as  special  instances,  corresponding  to  h = \ [Ju2  +v2  d£]2).  The  cubic  (5)  is 
the  simplest  (non-PH)  example  of  these  algebraically-rectifiable  curves,  with 
h(0  = (9fc2£2+4)3/2916fc4  — indeed,  it  is  the  unique  cubic  with  this  property. 
See  [53]  for  details  on  algebraically-rectifiable  quartics  and  quintics. 

Since  arc-length  parameterization  by  rational  functions  is  fundamentally 
impossible,  it  seems  natural  to  ask  “how  close”  we  can  approach  this  elusive 
ideal.  Consider,  for  example,  a degree-n  polynomial  curve  r(£):  a parameter 
transformation  £ € [ 0, 1 ] — * t € [ 0, 1 ] of  the  form 

t = {l-a)r 

a(l  — r)  + (1  — a )r 

gives  a rational  representation  of  the  same  degree,  and  offers  a single  degree  of 
freedom,  a,  to  control  the  “parameter  flow”  over  the  curve.  Using  the  integral 

1 = f (lr'(r)l  - l)2  dr 

Jo 


(7) 
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as  a measure  of  “closeness”  to  arc-length  parameterization  (for  which  1 = 0), 
the  value  of  a that  minimizes  (7)  can  be  found  [18]  as  the  unique  root  on  (0, 1) 
of  a quadratic  equation;  see  also  [36].  However,  this  optimal  parameterization 
offers  limited  scope  for  improvement,  since  we  fix  the  curve  degree  n. 

Another  approach,  based  on  the  polynomial  arc-length  functions  s(£)  of 
PH  curves,  employs  the  Legendre  series  to  compute  a convergent  sequence  of 
(constrained)  polynomial  approximations  £i(s),  ^(s), ...  to  the  inverse  of  this 
function  [19],  such  that 

lim  (k(s(0)  = 1 for  fe  [0,1], 

K— *00 

given  the  normalization  s £ [0, 1 ].  The  coefficients  of  (k{s)  can  be  determined 
through  closed-form  reduction  of  certain  integrals.  For  sufficiently  high  k , the 
re-parameterized  version  rjt(s)  = r(£/t(s))  comes  arbitrarily  close  to  the  exact 
arc-length  parameterization,  although  it  is  formally  of  degree  kn. 

§4.  Curves  from  Motion  I.  Centripetal  Forces 

By  integrating  the  description  of  the  forces  and  laws  of  motion  that  govern  a 
physical  system,  the  science  of  dynamics  provides  a rigorous  and  quantitative 
approach  to  analyzing  motions.  Perhaps  more  so  than  in  any  other  branch  of 
science  [12],  the  theoretical  canonization  and  empirical  triumph  of  dynamics 
are  the  fruits  of  a single  pre-eminent  mind:  Sir  Isaac  Newton. 

Perhaps  the  simplest  (non-trivial)  problem  of  dynamics  is  that  of  motion 
under  a centripetal  force  — i.e.,  a force  always  directed  toward  or  away  from 
a fixed  center,  whose  magnitude  depends  only  on  distance  r from  that  center. 
The  term  centripetal — “seeking  the  center”  — was  introduced  [60]  by  Newton 
in  his  De  motu  corporum  in  gyrum  of  1684  (in  recognition  of  the  fact  that,  to 
overcome  the  inertial  tendency  of  a body  to  move  in  a straight  line,  circular 
motion  requires  a steady  force  directed  toward  a fixed  center): 

I call  that,  by  which  a body  is  impelled  or  attracted  toward  some 
point  which  is  regarded  as  a center,  centripetal  force. 

The  basic  questions  concerning  centripetal  forces  are:  what  kinds  of  orbit  arise 
from  different  dependencies  of  the  force  on  r — and,  conversely,  knowing  the 
type  of  orbit,  can  we  deduce  the  dependence  of  the  force  on  r? 

In  the  late  17th  century  the  context  for  interest  in  such  questions  was,  of 
course,  the  search  for  an  explanation  of  Kepler’s  (empirical)  laws  of  planetary 
motion  — namely:  (i)  the  orbits  of  the  planets  are  ellipses,  with  the  sun  at  one 
focus;  (ii)  the  radial  line  between  the  sun  and  a planet  sweeps  out  area  at  a 
uniform  rate;  and  (iii)  the  squares  of  the  orbital  periods  are  proportional  to  the 
cubes  of  the  mean  distances  of  planets  from  the  sun.  As  we  now  know,  these 
are  direct  consequences  of  an  inverse-square  ( r~ 2)  gravitational  force.  Newton 
discovered  this  at  an  early  stage  in  his  career,  but  remained  characteristically 
secretive  about  it.  It  is  Edmond  Halley  who  deserves  credit  for  coaxing  Newton 
into  disseminating  his  arguments  and,  ultimately,  codifying  dynamics  through 
publication  of  the  Principia.  According  [13]  to  Abraham  de  Moivre: 
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In  1684  Dr  Halley  came  to  visit  him  at  Cambridge,  after  they 
had  been  some  time  together  the  Dr  asked  him  what  he  thought  the 
Curve  would  be  that  would  be  described  by  the  Planets  supposing 
the  force  of  attraction  towards  the  Sun  to  be  reciprocal  to  the  square 
of  their  distance  from  it.  Sr  Isaac  replied  immediately  it  would  be 
an  Ellipsis,  the  Dr  struck  with  joy  & amazement  asked  him  how  he 
knew  it,  why  saith  he,  I have  calculated  it,  whereupon  Dr  Halley 
asked  him  for  his  calculation  without  any  further  delay,  Sr  Isaac 
looked  among  his  papers  but  could  not  find  it,  but  he  promised  to 
renew  it,  & then  send  it  him  . . . 


In  fact  Halley,  as  Clerk  to  the  Royal  Society,  printed  the  Principia  at  his  own 
expense;  the  Society’s  funds  had  been  depleted  by  the  production  of  a Historia 
Piscium  ( History  of  Fishes)  that  failed  to  become  a best-seller.  Subsequently, 
Halley’s  salary  was  paid  entirely  in  copies  of  this  Historia  Piscium  [13]. 

The  inverse-square  nature  of  gravitational  force,  now  common  knowledge, 
was  established  by  Newton  as  the  rational  explanation  for  Kepler’s  laws.  It  is 
not  widely  known,  however,  that  the  Principia  thoroughly  analyzes  a variety  of 
power-law  (rn)  centripetal  forces,  and  shows  that  different  integer  exponents 
n yield  circular,  conic,  spiral,  and  other  orbits  in  an  often  surprising  manner. 
The  profundity  of  these  results  — and  their  anticipation  of  an  elegant  theory  of 
dual  centripetal  forces  due  to  K.  Bohlin  [7]  and  E.  Kasner  [38]  — has  recently 
been  emphasized  by  Arnol’d  [2],  Chandrasekhar  [10],  and  Needham  [45]. 

As  reflected  in  the  title  of  Chandrasekhar’s  recent  book  [10],  the  obscurity 
of  such  interesting  results  in  the  Principia  is  due  to  Newton’s  exclusive  reliance 
on  forbiddingly  Euclidean  argumentations.  The  modern  reader,  equipped  with 
predominantly  analytic/algebraic  skills,  is  usually  reduced  to  a state  of  dismay 
and  bewilderment  upon  a first  encounter  with  the  Principia.  Needham  [45] 
argues  convincingly  that  the  re-discovery  of  Newton’s  “geometrical  calculus” 
is  a very  rewarding  endeavor,  and  Arnol’d  [2]  gives  an  example  of  the  type  of 
problem  — the  transcendental  nature  of  certain  area  integrals  — that  seems 
“obvious”  to  Newtonian  thinking,  but  not  to  modern  modes  of  thought. 

To  elucidate  connections  between  the  Principia's  results  and  the  Bohlin- 
Kasner  theory  of  dual  centripetal  forces,  it  is  convenient  [2,44]  to  adopt  the 
complex-number  representation  z = x + i y = r e'e  of  the  Euclidean  plane. 
Under  a centripetal  force  proportional  to  the  n-th  power  of  distance  £from 
the  origin,  the  equation  of  motion  is  then 

Here,  the  centripetal  force  is  attractive  or  repulsive  according  to  whether  k is 
positive  or  negative.  Hooke’s  law,  for  example,  corresponds  to  n = +1,  while 
n = —2  represents  a Coulomb  (gravitational  or  electrostatic)  force. 

There  are  two  “constants  of  motion”  associated  with  solutions  z(t)  to  the 
differential  equation  (8)  — the  angular  momentum  and  energy, 


( dz\ 
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Fig.  3.  Constancy  of  angular  momentum  L under  an  r 2 centripetal  force. 

Although  we  may  formally  replace  jz|n+1/(n  + 1)  by  In \z\  when  n = —1,  this 
case  is  usually  excluded.  For  n < —2,  the  “potential  energy”  component  of  E 
is  naturally  negative,  and  tends  to  zero  as  |zj  — ► oo.  Conversely,  for  n > 0,  the 
potential  energy  is  naturally  positive  and  vanishes  when  \z\  — 0.  The  function 
ln|z|,  however,  diverges  for  both  \z\  —*  0 and  \z\  — * oo;  it  does  not  represent 
a satisfactory  scale-free  potential  energy  with  a natural  reference  value. 

The  polar  form  of  L given  in  (9)  serves  as  a reminder  of  its  geometrical 
interpretation,  namely,  the  rate  at  which  the  position  vector  sweeps  out  area  in 
the  orbit  — this  constancy  of  “areal  velocity”  is  expressed  by  Kepler’s  second 
law.  Section  II  of  Book  I of  the  Principia  is  devoted  to  “the  determination  of 
centripetal  forces,”  and  Newton’s  immediate  concern  [46]  is  to  prove: 

PROPOSITION  I.  THEOREM  I 

The  areas  which  revolving  bodies  describe  by  radii  drawn  to  an  immovable 
centre  of  force  do  lie  in  the  same  immovable  plane,  and  are  proportional  to 
the  times  in  which  they  are  described. 

PROPOSITION  II.  THEOREM  II 

Every  body  that  moves  in  any  curved  line  described  in  a plane,  and  by  a 
radius  drawn  to  a point  either  immovable,  or  moving  forwards  with  an  uniform 
rectilinear  motion,  describes  about  that  point  areas  proportional  to  the  times, 
is  urged  by  a centripetal  force  directed  to  that  point. 

In  other  words,  Newton  first  shows  that  “centripetal  force  L = constant” 
(this  principle  is  illustrated  in  Figure  3,  for  the  case  of  a Coulomb  r-2  force). 
Consider  now  the  conformal  map  z — > w of  the  plane  given  by 

w = za  , (10) 

under  which  the  orbit  z(t)  determined  by  equation  (8)  is  transformed  into  an 
orbit  tn(r),  where  t and  r denote  times  on  these  orbits  corresponding  to  an 
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angular  position  6 about  the  origins  of  the  z and  w planes.  Writing  z = re'e 
and  w = pe'e,  we  stipulate  that  these  orbits  have  equal  angular  momentum 

T o d6  o dd  . . 

L — r — = p — = constant , (11) 

dr  dr 

and  we  ask:  is  the  orbit  w(r)  also  the  solution  to  an  equation  of  motion 


d2w 

dr2 


(12) 


under  a power-law  centripetal  force  — and  if  so,  how  are  the  exponents  n,  m, 
a related?  It  transpires  that  this  problem  has  an  elegant  and  unique  solution: 
corresponding  force-law  exponents  n and  m are  related  by 


(n  + 3)(m  + 3)  = 4 , 


(13) 


and  the  exponent  a of  the  map  (10)  is  given  in  terms  of  them  by 

ti  + 3 2 

01  2 777  + 3 

To  derive  equations  (12)-(14),  we  first  note  from  (10)  and  (11)  that  derivatives 
with  respect  to  t and  r are  related  by 


iL  = |Z|2(1-“)1. 

dr  dt 


(15) 


Now  by  applying  (15)  to  to  = za  twice,  and  invoking  (8),  we  obtain 


d2717 

dr2 


2 a (1  — a) 
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dt 
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|z|n+1 

2(a-l) 


1 2— 3c* 


We  observe  that  by  choosing  2(a  — 1)  = 7i  + 1,  i.e.,  a = (n  + 3)/2  as  in  (14), 
the  expression  in  brackets  coincides  with  the  energy  constant  E in  (9).  With 
this  choice,  substitution  from  (10)  gives 

^ = - 2 a (a  - 1)  E [W|(2-3“)/“  . 

d Ti  |tl7| 

This  is  of  the  desired  form  (12),  with  k = 2a(a  — l)E  and  m = (2  — 3a)/a,  i.e., 
a = 2/(7n  + 3)  as  in  (14).  Finally,  equation  (13)  follows  from  the  individual 
relations  in  (14),  between  n and  a,  and  m and  a,  derived  above. 

Equation  (13)  describes  a hyperbolic  relation  between  the  force  exponents 
n and  m.  For  each  n,  except  —3,  the  orbit  z(t)  determined  by  (8)  is  mapped 
by  (10)  into  a dual  orbit  w(r)  of  equal  angular  momentum,  determined  by  (12). 
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Dual  orbits  corresponding  to  integer  force  exponents  are  of  special  interest  — 
they  are  (with  n > m)  as  follows: 


(a) 

n — +1, 

m = 

-2, 

z 

(b) 

T— ( 
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m = 
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T 

(d) 

n = —5, 

m = 

-5, 

z 

-+Z-1. 

Cases  (b)  and  (d)  identify  “self-dual”  forces  (as  noted  above,  however,  case  (b) 
is  usually  excluded  on  physical  grounds).  Case  (a)  reveals  the  beautiful  result 
that  an  orbit  z(t)  under  a linear  Hooke’s-law  force  is  mapped  by  tu  = z2  to 
an  orbit  w{r)  under  an  inverse-square  Coulomb  force.  Indeed,  Newton  shows 
that  the  Hooke  and  Coulomb  forms  are  the  only  centripetal  forces  that  admit 
conic  orbits.  In  this  regard,  see  [44]  for  an  interesting  anecdote  concerning  the 
£1  note  issued  to  commemorate  the  Principled s 300th  anniversary. 

To  investigate  the  geometry  of  orbits  under  power-law  centripetal  forces, 
it  is  convenient  to  employ  polar  coordinates  (r,  9).  Using  the  fact  that  r2#  = L, 
we  can  write  the  scalar  components  of  the  equation  of  motion  (8)  as 

r — L2r~3  + krn  = 0 and  2 r9  + r#  = 0 , (16) 


where  dots  denote  time  derivatives.  To  eliminate  the  time  variable,  and  obtain 
a purely  geometrical  description  of  the  orbit,  we  set  u = 1/r  and  note  that 


A = 0 -1  = Lu 
d t d 9 


2± 
de ' 


With  P = k/L2,  the  first  of  equations  (16)  can  then  be  transformed  [10,57]  to 
u"  + u - /3w"(n+2)  = 0,  (17) 


where  primes  denote  derivatives  with  respect  to  9.  By  solving  this  differential 
equation,  we  obtain  polar-coordinate  expressions  r(9)  = 1 /u(6)  describing  the 


shapes  of  orbits.  One  may  verify  the  functional  form  and  geometrical  nature 
of  some  representative  solutions  known  to  Newton: 

n = +1  : 

r(9)  = (/I  sin2  6 + cos2  0)-1/2 

ellipse,  center  at  origin; 

n = — 2 : 

r(9)  = (P  + cos#)-1 

ellipse,  focus  at  origin; 

n = —3  : 

r(9)  — exp  ( \JP  — 1 6) 

logarithmic  spiral; 

n = — 5 : 

r(6)  = \/PI‘2  cos  6 

circle  through  origin. 

The  cases  n = 1 and  -3  are  treated  in  PROPOSITION  X.  PROBLEM 
V and  PROPOSITION  IX.  PROBLEM  IV  of  the  Principia,  Book  I.  The  case 
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n = — 5 is  pathological,  since  the  orbit  passes  through  the  center  of  force!  In 
PROPOSITION  V.  PROBLEM  II,  Newton  actually  treats  a generalization  of 
the  n = — 5 case  — he  shows  that  a particle  p will  execute  a circular  orbit  if 
it  is  attracted  to  any  center  c by  a force  proportional  to  r~2£~3,  where  r is 
the  distance  of  p from  c,  and  t is  the  length  of  the  chord  containing  p and  c. 
When  c lies  on  the  circle,  we  have  £ = r,  and  hence  an  r~ 5 force. 

In  Section  III  of  Book  I,  Newton  is  concerned  with  “the  motion  of  bodies 
in  eccentric  conic  sections.”  He  treats  the  case  n = — 2 in  PROPOSITION  XI. 
PROBLEM  VI,  and  also  discusses  parabolic  and  hyperbolic  orbits.  Kepler’s 
third  law  is  also  derived,  in  PROPOSITION  XV.  THEOREM  VII. 

Incidentally,  the  case  n = 1 provides  the  correct  solution  to  the  problem 
of  motion  in  the  gravity  of  a permeable  rotating  Earth,  considered  by  Torricelli 
(see  §3).  If  the  Earth  is  a homogeneous  sphere  of  mass  M and  radius  FI,  the 
gravitational  force  at  distance  r from  the  center  is  equal  to 

GM{r/Rf 

r2 

G being  the  gravitational  constant  — i.e.,  it  is  proportional  to  r.  Thus,  the 
path  is  an  ellipse,  and  not  the  logarithmic  spiral  suggested  by  Torricelli  (which 
requires  an  r~3  force).  Actually,  it  is  a very  shallow  ellipse  — the  minor  axis  is 
smaller  than  the  major  axis  R by  the  dimensionless  factor  -J GM/oj2R 3 « 291, 
where  w = 2ir  rads/day  (this  factor  is  the  ratio  of  the  orbital  velocity  at  r = R 
to  the  tangential  velocity  uR  at  the  equator  due  to  the  Earth’s  rotation). 

A shallow  ellipse  is  an  obvious  perturbation  to  Galileo’s  simple  harmonic 
motion  through  a non-rotating  permeable  Earth,  and  in  retrospect  Torricelli’s 
conjectured  spiral  trajectory  — revived  by  Newton  in  1679  — may  seem  rather 
naive.  Newton  soon  redeemed  himself,  however,  through  the  publication  of  his 
Principia  in  1686,  which  contains  the  correct  solution  as  part  of  a remarkably 
comprehensive  theory  of  orbital  motions  under  centripetal  forces. 

It  is  a sobering  experience,  for  the  modern  reader,  to  pierce  the  Principia' s 
veil  of  geometrical  argumentations,  and  appreciate  its  profound  insights.  To 
contemporaries,  Newton’s  creation  was  a virtually  miraculous  event  — Halley 
composed  an  ode  to  preface  the  Principia,  extolling  “the  illustrious  man”  and 
his  work,  “a  signal  distinction  of  our  time  and  race”  [46]: 

Matters  that  vexed  the  minds  of  ancient  seers, 

And  for  our  learned  doctors  often  led 
To  loud  and  vain  contention,  now  are  seen 
In  reason ’s  light,  the  clouds  of  ignorance 
Dispelled  at  last  by  science  . . . 

However,  Newton  was  not  exempt  from  the  sarcasm  of  critics,  such  as  the  poet 
Alexander  Pope  [51],  who  were  often  more  eloquent  in  their  converse  views: 

Superior  beings,  when  of  late  they  saw 
A mortal  man  unfold  all  Nature’s  law, 

Admired  such  wisdom  in  an  Earthly  shape, 

And  showed  a NEWTON  as  we  show  an  ape. 
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To  conclude,  we  note  that  the  singular  orbit  r(0)  = \//3/2  cos  6,  identified 
above  for  an  r~5  force,  corresponds  to  the  case  where  the  orbital  energy 

E ~ 2 + 2^  ~ 4^ 

is  zero:  it  is  the  analog  of  zero-energy  parabolic  orbits  under  an  r~2  force  (see 
Figure  4).  As  r — ► 0 (t  — > nk/8L3),  the  positive  kinetic  energy  and  negative 
potential  energy  both  become  infinite  in  a manner  such  as  to  maintain  E — 0. 
The  angular  momentum  L is  also  conserved,  in  a limiting  sense,  as  r -+  0. 
For  orbits  with  E / 0,  equation  (17)  can  be  integrated  to  obtain 


where  7 = E/L2  and  r = r0  for  # = 0.  A further  reduction,  giving  r explicitly 
in  terms  of  6,  is  possible  upon  introducing  Jacobian  elliptic  functions  [41],  but 
we  shall  not  pursue  it  here.  A special  case  is  a circular  orbit  u = uo  = l/r0,  of 
energy  E = \(2L2  — ku{ \)u2,  but  this  is  highly  unstable  - — any  perturbation 
will  cause  r to  rapidly  decay  to  0 or  grow  to  00  (see  Figure  5,  comparing  orbits 
under  r-5  and  r~ 2 forces  perturbed  by  introducing  an  initial  negative/positive 
radial  velocity  corresponding  to  a 10-3  fractional  change  in  E). 
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In  fact,  only  the  n = +1  and  —2  force  laws  admit  stable  periodic  orbits, 
in  the  sense  that  perturbations  to  them  always  produce  similar,  “neighboring” 
closed  orbits  (a  proof  is  given  in  Appendix  A of  Goldstein  [28]). 

§5.  Curves  from  Motion  II.  Four-bar  Linkage 

A mechanism  is  a device  that  is  designed  to  transform  “input”  motions  and 
forces,  from  a given  power  source,  into  “output”  motions  and  forces  — better 
suited  for  use  in  some  practical  application.  A mechanism  typically  comprises 
several  rigid  members  connected  by  joints  that  allow  certain  types  of  relative 
motion.  According  to  the  Kempe  theorem  [39],  mechanisms  that  employ  only 
revolute  and  prismatic  joints  can  (in  principle)  be  designed  to  produce  motion 
along  any  plane  algebraic  curve.  We  consider  here  the  coupler  curves  of  planar 
four-bar  linkages,  which  serve  to  transform  a rotational  input  motion  into  an 
output  motion  along  some  general  curved  trajectory. 

(In  early  studies,  such  mechanisms  were  called  three- bar  linkages,  since 
the  “ground”  link  was  not  counted.  It  is  now  customary  to  include  it,  to  give 
a closed  kinematic  chain.  The  idea  of  kinematic  chains  was  introduced  by  the 
German  engineer  Franz  Reuleaux,  in  his  Theoretische  Kinematik  of  1875). 

Figure  6 shows  the  geometrical  configuration  of  a four-bar  linkage.  Such 
mechanisms  are  found  in  diverse  contexts  (windshield  wipers,  electric  shavers, 
cranes,  etc.).  Historically,  the  most  famous  example  was  the  “parallel  motion” 
mechanism  devised  by  James  Watt  (1736-1819)  for  his  double-acting  steam 
engine  of  1782.  In  Figure  6,  the  link  of  length  k is  held  fixed,  while  links  of 
length  r and  R pivot  about  its  two  endpoints.  These  links  are  connected  by 
a further  link  of  length  c,  whose  ends  are  thus  constrained  to  lie  on  circles  of 
radii  r and  R centered  on  (0, 0)  and  (k,  0).  A point  at  a fixed  position  relative 
to  the  link  of  length  c thus  traces  a locus,  called  the  coupler  curve,  when  the 
links  of  length  r and  R rotate  (see  Figure  7).  We  identify  a specific  point  p 
by  taking  it  to  be  the  apex  of  a triangle  of  sides  a and  b,  with  base  c. 
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Fig.  7.  Generation  of  a coupler  curve  by  a four-bar  linkage. 


According  to  the  Grashof  theorem  [29],  the  sum  of  lengths  of  the  shortest 
and  longest  links  should  not  exceed  the  sum  of  lengths  of  the  other  two  links, 
if  there  is  to  be  continuous  relative  rotation  between  two  links.  This  condition 
is  satisfied  by  the  configuration  shown  in  Figure  6,  a crank-rocker  mechanism: 
the  link  of  length  r (the  crank)  completes  full  revolutions,  while  that  of  length 
R (the  rocker)  oscillates  through  partial  revolutions;  the  link  of  length  c is  the 
coupler.  Other  possible  mechanisms  — the  double-crank  or  double-rocker  — 
are  obtained  by  varying  the  lengths  k,  r,  R,  and  c of  the  links. 

Four-bar  linkages  are  capable  of  generating  a rich  variety  of  curved  paths: 
the  Hrones-Nelson  “atlas”  [32],  for  example,  illustrates  over  seven  thousand 
different  forms  of  the  coupler  curve!  These  paths  are  all  (parts  of)  an  algebraic 
curve  of  degree  6,  that  depends  on  six  parameters  — the  dimensions  k,  r,  R, 
a,  b,  and  c.  Its  equation  can  be  succinctly  expressed  [4]  in  the  form 

f(x,y)  = u2(x,y)  + v2(x,y ) - w2(x,y)  = 0 (18) 


where 

u(x,  y)  — a [ (x  — k)  cos  7 + y sin  7]  (x2  + y2  + b2  — r2) 

— bx  [ (x  — k )2  + y2  + a2  - R2  ] , 

v(x,y)  = a [ (x  — k)  sin7  — 7/COS7]  (x2  + y2  + b2  — r2)  (19) 

+ by  [(x  - k)2  + y2  + a2  - R2 } , 
w(x,y)  = 2ab  [ x(x  - k ) sin  7 + y2  sin  7 — ky  cos  7 ] , 

and  we  have  introduced  the  angle  7 = cos-1  (a2  + b2  — c2)/2ab  in  lieu  of  c.  The 
curve  defined  by  (18)  and  (19)  has  an  ordinary  triple  point  at  each  of  the  two 
circular  points  at  infinity,  and  three  affine  double  points  (two  of  which  may 
be  complex  conjugates)  that  always  lie  [4]  on  the  circle 

x2  — kx  + y2  — ky  cot  7 = 0, 


(20) 
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illustrated  in  the  example  in  Figure  9 below.  The  coupler  curve  falls  one  short 
of  the  maximum  of  10  double  points  that  an  algebraic  curve  of  degree  6 may 
have,  and  is  therefore  of  genus  1 — i.e.,  it  is  an  elliptic  curve. 

For  a crank-rocker,  the  curve  defined  by  (18)  and  (19)  always  comprises 
two  real  loops.  The  physical  mechanism  traces  just  one  of  them:  to  trace  the 
other  loop,  the  initial  configuration  of  the  linkage  must  be  changed.  Further 
examples  of  crank-rocker  coupler  curves  are  shown  in  Figure  8.  For  double- 
crank and  double-rocker  mechanisms,  the  coupler  curve  also  has  two  loops,  but 
in  the  latter  case  the  mechanism  cannot  trace  either  loop  entirely.  Four-bar 
linkages  that  do  not  satisfy  the  Grashof  condition  exhibit  single-loop  coupler 
curves,  which  may  self-intersect  (as  in  the  Figure  9 example  — note  that  the 
mechanism  cannot  trace  the  entire  curve).  Equations  (18)  and  (19)  encompass 
all  these  forms  for  suitable  choices  of  k,  r,  R,  a,  b,  and  c. 


Fig.  9.  Single-loop  coupler  curve  for  a non-Grashof  mechanism. 
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Fig.  10.  Peaucellier  mechanism:  circular  motion  of  P yields  linear  motion  of  P' . 


A remarkable  property  of  four-bar  linkage  coupler  curves  is  expressed  by 
the  Roberts-Chebyshev  theorem  [50]  — a given  coupler  curve  may  actually  be 
traced  by  three  different  four-bar  mechanisms  (which  are  said  to  be  cognates 
with  respect  to  the  given  curve).  This  is  not  merely  a mathematical  curiosity: 
once  a path  is  realized  through  a specific  mechanism,  one  of  its  cognates  may 
be  found  to  produce  better  force  transmission  characteristics. 

The  design  of  a mechanism  for  a given  path  may  be  based  on  consultation 
of  an  “atlas”  of  coupler  curves  [32],  or  use  of  numerical  methods  to  find  the 
linkage  dimensions  that  will  give  a locus  interpolating  prescribed  points  [6,58]. 
Allowing  for  freedoms  in  the  choice  of  coordinate  system,  the  general  coupler 
curve  defined  by  (18)  and  (19)  can  be  made  to  interpolate  nine  points,  but 
determining  the  mechanism  parameters  involves  solving  a formidable  system  of 
non-linear  algebraic  equations.  Wampler  et  al.  [58]  have  shown,  for  example, 
that  (counting  cognates)  the  nine-point  problem  has  4326  solutions,  many  of 
which  may  correspond  to  complex  values  for  the  link  dimensions,  or  interpolate 
the  discrete  points  on  incompatible  portions  of  the  coupler  curve. 

As  previously  noted,  one  of  the  first  applications  of  four-bar  linkages  was 
to  transform  the  reciprocating  linear  motion  of  a piston  into  rotary  motion  of  a 
shaft.  Watt,  Chebyshev,  Roberts,  and  others  proposed  approximate  solutions 
to  this  problem,  but  linkages  that  offer  exact  transformations  between  linear 
and  circular  motion  were  not  known  until  1864,  when  a captain  in  the  French 
army,  A.  Peaucellier,  devised  the  mechanism  shown  in  Figure  10,  comprising 
four  links  of  length  a and  two  of  length  b (>  a).  If  the  point  P is  constrained 
to  move  on  a circle  passing  through  the  pivot  0,  the  point  P'  traces  a straight 
line,  such  that  O,  P , P'  always  remain  collinear,  and  the  distances  r and  r' 
of  P and  P'  from  O satisfy  rr'  = b2  — a2.  Thus,  P and  P'  are  images  of  each 
other  under  inversion  in  the  circle  with  center  0 and  radius  y/b 2 — a2. 

Further  details  on  the  geometrical  and  kinematical  properties  of  coupler 
curves  may  be  found  in  standard  texts  [4,16,30,33]  on  kinematics  — see  also 
[47]  for  an  interesting  history  of  coupler  curve  synthesis  methods. 
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§6.  Motion  from  Curves  I.  Multi-axis  CNC  Machines 

We  have  described  above  how  interesting  curves  can  arise  from  motions  under 
specified  kinematical  or  dynamical  constraints.  Computer  numerical  control 
(CNC)  technology  is  concerned  with  the  converse  problem  — i.e.,  the  physical 
realization  (by  cutting  tools,  robot  arms,  sensors,  etc.)  of  motions  specified 
by  geometrical  paths  and  given  speeds  or  feedrates  along  these  paths. 

To  produce  a desired  motion,  a CNC  machine  must  drive  each  of  its  axes 
in  an  independent  but  coordinated  manner.  The  controller  algorithm  employs 
digital  representations  of  space  and  time:  the  unit  of  time,  or  sampling  interval 
(typically  rsj  1 millisecond),  is  defined  by  a “clock”  within  the  algorithm,  while 
the  basic  length  unit  (BLU,  typically  ~ 10  microns),  or  spatial  resolution  of 
the  machine,  is  determined  by  position  encoders  mounted  on  its  axes. 

Within  each  sampling  interval  At,  the  controller  must  compare  the  actual 
position  of  each  axis  (as  measured  by  the  encoders)  with  the  intended  position 
(computed  from  the  specified  paths  and  feedrates  by  a real-time  interpolator). 
The  discrepancy  between  the  actual  and  desired  positions  is  used  to  generate 
control  signals  for  the  machine  drives,  ensuring  that  the  specified  paths/speeds 
are  accurately  realized.  The  discrete  positions  on  a curved  path  r(£)  computed 
by  the  interpolator  are  known  as  reference  points  — they  are  identified  by  the 
sequence  £oi  £ii  62  > ■ • • of  parameter  values  satisfying  £0  = 0 and 

^ = At  for  * = 1,2,...,  (21) 

J(k-i  K 


where  a(f)  = |r'(£)|  is  the  parametric  speed  of  the  curve  and  V is  the  (constant 
or  variable)  feedrate.  Since  the  integral  does  not  ordinarily  have  a closed-form 
reduction,  and  the  unknowns  are  limits  of  integration,  equation  (21)  is  difficult 
to  solve  accurately  and  efficiently  — even  if  V is  constant. 

Because  of  this  computational  difficulty,  it  is  customary  to  employ  simple 
(piecewise-linear/circular)  “G  code”  approximations  to  curved  tool  paths  [1]. 
Compared  to  its  electromechanical  hardware  sophistication,  the  part  program 
data  that  drives  a CNC  machine  is  embarrassingly  crude.  Some  authors  [11,61] 
have  proposed  to  drive  CNC  machines  along  general  curved  paths  by  invoking 
approximate  solutions  to  (21),  based  on  the  Taylor-series  expansion 


& = 6> 


V 

-1  + — A t + 


Kv'-r^v) 


(At)2 


+ 


(22) 


where  primes  indicate  derivatives  with  respect  to  the  curve  parameter  £,  and 
it  is  understood  that  a,  r',  r",  V,  V',  etc.,  are  evaluated  at  £k-i- 

The  extension  of  (22)  to  cubic  and  higher-order  terms  incurs  complicated 
coefficients,  and  is  thus  ill-suited  to  real-time  computation.  Truncation  errors 
are  inevitable  with  this  approach  (most  implementations,  in  fact,  retain  only 
the  linear  term).  Note  also  that,  for  a non-constant  feedrate,  the  variation  of 
V cannot  be  usefully  specified  as  a function  of  £.  It  must  be  given  in  terms  of 
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Fig.  11.  Curvature  dependent  feedrate  (24)  for  a constant  material  removal  rate. 


a physically  meaningful  variable,  such  as  time  t,  arc  length  s,  or  curvature  k 
— in  terms  these  variables,  the  derivative  of  V in  (22)  may  be  written  as 

V'  - - — - — - — — 

V df  ds  ds  d/t 

The  Pythagorean-hodograph  curves  offer  an  elegant  and  rigorous  solution 
to  this  dilemma  [24],  For  PH  curves,  the  integral  (21)  admits  a simple  analytic 
reduction  — not  only  for  constant  feedrate,  but  also  varying  feedrates  specified 
in  a number  of  useful  ways,  such  as: 

(a)  any  function  V(t)  of  time  with  known  indefinite  integral; 

(b)  a linear  or  quadratic  polynomial  F(s)  in  the  arc  length ; 

(c)  simple  rational  expressions  V (k)  in  the  local  curvature ; 

(d)  constant  feedrate  V along  an  offset  to  a specified  curve. 

In  (almost)  all  these  cases,  the  interpolation  equation  (21)  reduces  to  the  form 

«(&)  = *’(•..,&- 1),  (23) 


where  s(£)  is  the  polynomial  arc-length  function  of  the  PH  curve,  and  F is  a 
known  elementary  function  of  the  parameters  describing  the  feedrate  variation 
and  the  preceding  reference  point  £k-i-  Since  s(£)  is  a monotone  polynomial, 
equation  (23)  has  a unique  real  root  for  the  value  of  £*,,  which  may  be  obtained 
to  machine  precision  by  a few  Newton-Raphson  iterations  starting  from 

The  ability  to  perform  real-time  interpolation  with  continuously  varying 
feedrates  is  extremely  useful  in  a variety  of  practical  problems.  For  example, 
the  curvature-dependent  feedrate 


1 + K,(d  — |6) 


(24) 
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time  (teconds)  length  (mm) 

Fig.  12.  Measured  feedrate  (left)  and  cutting  force  (right)  using  the  function  (24). 

can  suppress  machining  force  variations  [20]  when  a fixed  depth  of  cut  6 is  to 
be  removed  along  a curved  path  r(£)  with  a tool  of  radius  d (see  Figure  11:  the 
material  removal  rate  is  higher  at  the  “concave”  location  b than  the  “convex” 
location  c if  a constant  feedrate  Vo  is  employed).  Figure  12  shows  the  measured 
feedrate  and  (time-averaged)  cutting  force,  as  obtained  from  an  experimental 
implementation  of  the  PH  curve  interpolator  for  (24)  on  a CNC  mill. 

An  important  use  of  time-dependent  feedrates  V (t)  is  the  specification  of 
smooth  accelerations  and  decelerations  along  curved  paths.  This  is  especially 
important  for  high-speed  machining,  in  which  the  dynamical  issues  of  starting 
and  stopping  high-speed  motions  on  curved  paths  become  a serious  concern. 
With  G code  part  programs,  for  example,  acceleration/deceleration  intervals 
may  span  many  short  linear/circular  segments,  and  thus  require  a cumbersome 
real-time  block  “look  ahead”  capability  for  their  implementation. 

Using  PH  curves,  on  the  other  hand,  it  is  easy  to  specify  a smooth  feed 
acceleration  from  feedrate  V = 0 for  t < 0 up  to  a desired  constant  feedrate 
Vm  for  t > T along  a curve  r(£).  If  r = t/T  is  the  “normalized”  time  during 
the  acceleration  interval  t € [0,  T],  we  use  the  polynomial  feedrate  function 

^)  = EF‘Q(i-r)n'M  (25) 

fc=0  ' ' 

of  odd  degree  n,  with  V0  = • • • = Vr(n_1)/2  = 0 and  V(„+1)/2  = • • • = Vn  = Vm. 
This  gives  Ct'n~1^2  continuity  with  V = 0 for  t < 0 and  V = Vm  for  t > T — 
in  particular,  n = 3 and  n = 5 yield  C1  and  C2  feedrate  variations. 

The  interpolation  equation  for  such  time-dependent  feedrate  functions 
on  PH  curves  is  remarkably  simple  — it  is  precisely  of  the  form  (23),  with  the 
right-hand  side  being  simply  the  integral  F(t)  of  (25),  a polynomial  of  degree 
n + 1.  Only  the  right-hand  side  “constant”  changes  on  using  higher  degrees  n, 
and  the  incremental  cost  of  evaluating  this  constant,  in  each  sampling  interval 
At,  is  insignificant  for  any  “reasonable”  choice  of  the  degree  n. 

There  are  many  other  possibilites  for  specifying  and  optimizing  feedrates 
along  PH  curves.  For  example,  curved  paths  are  realized  on  CNC  machines  by 
coordinated  motions  of  independently-powered  axes,  and  the  chosen  feedrates 
and  feed  accelerations  should  not  impose  demands  on  the  motors  of  each  axis 
that  exceed  their  torque  or  power  capacity.  For  PH  curves,  a thorough  analysis 
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of  this  problem  is  possible  [25],  allowing  an  a priori  determination  of  safe  fixed 
feedrates  and  minimum  feed  acceleration  intervals  for  a given  path  geometry. 


§7.  Motion  from  Curves  II.  Quaternion  Methods 

We  have  only  been  concerned,  thus  far,  with  the  motion  of  points  in  Euclidean 
spaces.  The  problems  of  motion  in  non-Euclidean  spaces,  or  motion  of  bodies 
of  finite  extent  (involving  changes  of  position  and  orientation),  are  much  more 
challenging  — they  arise  frequently  in  animation,  robotics,  5-axis  machining, 
dynamics,  and  many  other  applications.  Thus,  to  conclude,  we  briefly  consider 
some  basic  problems  in  the  use  of  quaternions  [5,8,52]  to  formulate  spatial  rigid 
body  motions  as  time-parameterized  loci  in  non-Euclidean  spaces. 

This  subject  has  seen  intense  interest  [17,27,35,37,40,48,54,59]  — see  also 
the  extensive  bibliography  in  [52]  — in  recent  years,  and  substantial  progress  in 
“motion  design”  has  been  made.  Nevertheless,  the  fundamental  problems  we 
encounter  with  purely  translational  (point)  motion  carry  over  to  and,  indeed, 
become  much  deeper  in  the  context  of  general  spatial  motions. 

The  revival  of  quaternions  in  motion-design  applications  helps  remedy  a 
steady  historical  decline  of  interest  in  them.  The  introduction  of  quaternions 
by  Sir  William  Hamilton  (1805-1865)  predates  (and  subsumes)  development 
of  the  “ordinary”  vector  analysis  in  1R3  by  James  Clerk  Maxwell  (1831-1879), 
Josiah  Willard  Gibbs  (1839-1903),  and  Oliver  Heaviside  (1850-1925),  who  — 
along  with  later  generations  of  physicists  — considered  quaternions  to  be  an 
unduly  cumbersome  medium  for  describing  the  laws  of  nature  [14]. 

A general  displacement  in  Euclidean  3-space  H3  can  be  interpreted  as  a 
screw  displacement  — i.e.,  a rotation  about  a fixed  axis  and  a translation  along 
that  axis  [8].  Six  parameters  are  required  to  describe  such  displacements.  Let 
P — ( X , Y,  Z)T  and  p = ( x , y,  z)T  be  point  coordinates  in  a “fixed”  frame  E 
and  a “movable”  frame  a.  The  spatial  displacement  carrying  E into  a may 
be  described  by  an  orthogonal  rotation  matrix  M and  a translation  vector  d: 

P = Mp+  d. 


The  rotation  matrix  can  be  expressed  in  terms  of  Euler  parameters 


co  = cos  \(j> , cj  = Asini^>,  C2  = p sin  \<j> , C3  = ^sin^, 
where  (A,  p,  v)  are  direction  cosines  of  the  axis  and  </>  is  the  rotation  angle,  as 


M = 


r2  1 r2  2 _ 2 

c0  ' C1  c2  c3 

2(c2Ci  + C0C3) 
2(C3C!  — C0C2) 


2(ciC2  - c0c3) 

r2  _ .2  1 .2  _ „2 

c0  C1  ' c2  c3 
2(c3C2  + C0Ci) 


2(ciC3  + coc2) 
2(c2c3  - c0ci) 

r2  r2  r2  1 _2 

c0  C1  c2  T C3 


Note  that,  since  the  Euler  parameters  satisfy  the  normalization  condition 

4 + ci  + 4 + c3  = 1 > (26) 
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only  three  are  independent.  In  lieu  of  the  translation  vector  d = (di , d2.,d3), 
Study  [56]  introduced  four  new  parameters: 

c'o  = (cidi  + C2CZ2  + c3d3)/2 , 
dx  = (— codi  — 03^2  + c2d3)/2 , 
c2  = (c3^1  — co^2  — c\d3)/2 , 

C3  = {—c2d\  -f  cid2  — cod3)/2 , 

which,  by  definition,  satisfy  the  constraint 


coCq  + £\c'i  + C2C2  + C3C3  = 0 . (27) 

The  set  of  all  displacements  in  1R3  can  then  be  regarded  as  points,  or  soma,  in 
a 6-dimensional  space  spanned  by  eight  coordinates  (c0,  ci,  C2,  C3,  Cq,  c\ , c'2 , C3) 
subject  to  the  two  algebraic  constraints  (26)  and  (27). 

A compact  and  elegant  algebraic  description  of  these  soma  is  obtained  [8] 
by  combining  the  eight  coordinates  into  a dual  quaternion  of  the  form 

C = Co  + ecg  + i (ci  + eci)  + j (02  + ec^)  + k (03  + ec'3)  . (28) 

Here,  the  quaternion  basis  elements  satisfy  the  multiplication  rules 

ij  = k,  jk  = i,  ki=j,  i2=j2  = k2  = -l 

(so  that  multiplication  is  non-commutative:  ji  = — ij,  etc).  The  components 
of  the  quaternion  (28)  are  “dual  numbers”  of  the  form  x + ex'  for  real  x,  x' 
— the  dual  basis  element  e satisfies  e2  = 0 (^  e).  The  relations  (26),  (27)  and 
E2  = 0 ensure  that  (28)  is  a unit  dual  quaternion:  its  components  satisfy 

(co  + £Cg)2  + (ci+ec^)2  + (02  + ec^)2  + (c3  + ec'3)2  = 1. 

For  a spatial  displacement  specified  by  the  unit  dual  quaternion  (28),  we  can 
extract  the  geometrical  parameters  as  follows: 

d = 2 (c[,c  — coc'  — c x c') , <f>  = 2 cos-1  co  , (A,  fi,  v)  = — — , 

sin 

where  c = (cj,  02,03),  c'  = (c'^c^c'g),  0 < cos-1  co  < 7r,  and  x is  the  familiar 
vector  cross  product  in  1R3.  Note  also  that  the  translation  distance  is 

d = |d|  = 2-\Jc'q  + c'j2  + Cj2  + C32  . 

When  eg  = c(  = c'2  = C3  = 0 the  dual  quaternion  (28)  specifies  a pure  rotation 
(the  non-commutativity  of  quaternion  multiplication  reflects  the  importance 
of  the  order  in  which  spatial  rotations  are  executed  — see  Figure  13). 
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Fig.  13.  Effect  of  rotations  Rx(a)Rz(f3)  and  Rz(0)Rx(oi)  applied  to  a vector  v. 

The  principal  advantage  of  the  form  (28)  is  that  the  outcome  of  successive 
displacements,  A followed  by  B,  corresponds  to  the  ordered  product  C = BA 
of  their  dual  quaternion  representations.  With 

A = ag  -+•  ea'0  + i (aj  + sa'i)  + j (a?  + £02)  d k (°3  d £<*3) , 

B = b0  + eb'0  + i (61  + £b[)  + j (62  d £62)  + k (63  + £63) , 

the  elements  of  C are  homogeneous  quadratic  forms  in  those  of  A and  B: 


Co  = ao&o  — aibi  — O2&2  — 0363  ! 

ci  = aibo  + aobi  + CI3&2  — 0263  , 

C2  = 0260  + Oo&2  + «1&3  — a3&l  > 

C3  = 0360  + d0^3  + 02&1  — a1^2  i 

Cq  = a06o  + a'obo  — a.\b'i  — a'lbi  — — a'2b 2 — 0363  — (1363  , 

Cj  = ai&Q  + a\bo  + aob[  + a'gb  1 + 0363  + 0362  — 02^3  — ®2^3  j 

c2  026(5  d a(>6o  d Go ^2  d 0(562  d G163  + 0^63  — 036^  — O 3 6 1 , 

C3  — 036(5  d G360  d 0(563  + 0(563  + 02^1  + O261  — 01&2  — O162  • 

Now  a general  rigid  body  motion,  involving  both  translation  and  rotation, 
corresponds  to  a locus  of  points  in  soma  space,  where  a value  of  the  time  t is 
associated  with  each  point.  It  is  tempting  to  invoke  a unit  dual  quaternion 
C(t)  parameterized  explicitly  by  “simple”  — i.e.,  (piecewise)  polynomial  or 
rational  — functions  of  time,  and  most  motion  design  schemes  rely  upon  this 
model.  Typically,  a sequence  Ci, . . . , C/v  of  displacements  at  times  t\, ...  ,tN 
are  specified,  and  one  seeks  a “smooth”  motion  interpolating  them. 

In  §2  we  emphasized  that  motion  specification  is  as  much  concerned  with 
velocities  and  accelerations  — determined  by  the  precise  nature  of  the  time 
parameterization  — as  with  the  path  geometry,  and  this  precept  also  holds  for 
spatial  rigid  body  motions.  An  ad  hoc  or  “indirect”  time  parameterization  of  a 
quaternion  locus  C(t)  can  incur  linear/angular  velocities  or  accelerations  that 
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are  undesirable  or,  at  least,  only  determinable  a posteriori  — after  the  motion 
has  been  specified,  rather  than  being  an  integral  part  of  its  specification. 

For  rational  point  motion,  the  difficulty  in  simultaneously  specifying  both 
path  geometry  and  speed  along  the  path  arises  from  the  fact  that  curves  do 
not,  in  general,  admit  rational  arc  length  representations.  As  indicated  in  §6, 
the  use  of  special  (PH)  curves  can  resolve  this  problem.  For  rational  rigid  body 
motions,  however,  the  difficulty  is  not  just  one  of  computation,  but  also  the 
more  fundamental  issue  of  how  we  characterize  “distance  travelled”  in  terms  of 
both  the  translational  and  orientational  components.  Should  they  be  treated 
together,  or  separately?  In  other  words,  can  we  introduce  a suitable  metric  for 
soma  space  that  allows  us  to  define  “arc  length”  along  a unit  dual  quaternion 
locus  C(f),  and  thus  formulate  methods  to  specify  both  positions/orientations 
and  linear /angular  speeds  for  rational  spatial  motions? 

Ravani  and  Roth  [49]  have  proposed,  by  analogy  with  elliptic  geometry,  a 
metric  that  yields  a dual  number  value  for  distances  in  soma  space.  However, 
the  use  of  this  metric  in  motion  design,  or  of  alternate  (real-valued)  functions 
that  exhibit  the  usual  properties  of  metrics,  remains  to  be  explored. 

§8.  Closure 

With  its  opening  AXIOMS,  or  LAWS  OF  MOTION,  the  Principia  establishes 
uniform  motion  as  the  natural  state  of  a free  body.  The  forces  that  act  upon 
bodies  incur  deviations  from  uniform  motion,  in  a deterministic  though  subtle 
manner  that  reveals  appealing  and  useful  connections  between  geometry  and 
kinematics.  Although,  in  this  rather  brief  and  eclectic  survey,  we  have  offered 
only  a few  illustrative  anecdotes  on  this  theme,  we  hope  they  have  stirred  the 
interest  of  some  inquisitive  readers,  and  have  thus  helped  to  promote  further 
theoretical  developments  and  practical  applications. 
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